Extreme black hole simulations: 
collisions of unequal mass black holes and the point particle limit 
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Numerical relativity has seen incredible progress in the last years, and is being applied with success 
to a variety of physical phenomena, from gravitational-wave research and relativistic astrophysics 
to cosmology and high-energy physics. Here we probe the limits of current numerical setups, by 
studying collisions of unequal mass, non-rotating black holes of mass-ratios up to 1:100 and making 
contact with a classical calculation in General Relativity: the infall of a point-like particle into a 
massive black hole. 

Our results agree well with the predictions coming from linearized calculations of the infall of 
point-like particles into non-rotating black holes. In particular, in the limit that one hole is much 
smaller than the other, and the infall starts from an infinite initial separation, we recover the point- 
particle limit. Thus, numerical relativity is able to bridge the gap between fully non-linear dynamics 
and linearized approximations, which may have important applications. Finally, we also comment 
on the "spurious" radiation content in the initial data and the linearized predictions. 

PACS numbers: 04.25.D-, 04.25. dg, 04.25.Nx, 04.30.-w 
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I. INTRODUCTION 

In recent decades, black holes have started playing a 
key role in a variety of processes in astrophysics, gravita- 
tional wave physics and high-energy physics. Following 
the 2005 breakthroughs ll|-y|: numerical relativity has 
been an essential tool in the modeling of black-hole bina- 
ries in the strong-field regime. At the same time it has 
become clear that detailed studies of black-hole systems 
often involve a close interplay between fully non-linear 
numerical simulations and approximation techniques of 
various types. For example, the generation of gravita- 
tional wave (GW) template banks for use in the analy- 
sis of observational data from laser intcrferometric GW 
detectors LIGO, VIRGO, GEO600, LCGT or LISA re- 
quires the combination of numerical relativity with post- 
Newtonian or other techniques; sec Refs. p|-[lO| and ref- 
erences therein. Post- Newtonian studies have also played 
an important role in the guidance of the numerical in- 
vestigation of the black-hole recoil, most notably in the 
discovery of the so-called superkicks and its possible sup- 
pression due to spin alignment In the context of 
high-energy collisions of black holes, linearization tools 
such as the zero- frequency limit or point-particle calcula- 
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tions provide valuable insight into the scattering thresh- 
old and GW emission of black- hole collisions in four and 
higher-dimensional spacetimes [l6j . A particular class of 
black-hole binaries of high relevance for the spaceborne 
LISA (or a similar future spaceborne) observatory, the 
so-called extreme-mass-ratio inspirals, represent a par- 
ticularly difficult challenge to numerical relativity and 
their modeling relies heavily on perturbative methods 
and self-force calculations; see Refs. [T7l - l22j and refer- 
ences therein. 

With the above as motivation, it is vital to obtain a 
detailed understanding of the range of validity of the var- 
ious types of approximation methods. At the same time, 
these methods provide valuable tools to calibrate the ac- 
curacy of numerically generated solutions to the Einstein 
equations. The purpose of this paper is to provide such 
a study for the case of a classical calculation in general 
relativity, the head-on infall of a point-particle (PP) into 
a black hole [H. 

In recent years, numerical relativity has started prob- 
ing the intermediate mass-ratio regime by evolving the 
final orbits of (approximately) quasi-circular inspirals of 
black-hole binaries with mass-ratio q = m2/mi = 1/10 
[U, [2^ ; by comparing numerical results with perturba- 
tive calculations employing the fully numerical black-hole 
trajectories for mass ratios up to g = 1/20 (2^ and most 
recently, the first numerical evolution of a black-hole bi- 
nary with q = 1/100 (27j . In this work, we restrict our 
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attention on the head-on hmit of the colhsion of black 
holes, for two reasons: (i) the lower computational cost 
due to the higher degree of spacetime symmetry and the 
absence of the lengthy inspiral phase and (ii) the avail- 
ability of high-precision results in the PP limit. 

In our study we will make extensive use of the calcu- 
lation by Davis et al. [1^ who model in the PP limit the 
collision of a small object of mass m with a black hole 
of mass M ^ m. In the original calculation the particle 
was falling from rest at infinity, and the total radiated 
energy was found to be 
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This setting has been generalized to arbitrary initial dis- 
tance and boost, in which case initial data and conse- 
quent spurious radiation play a role [H, [28l - [3H . 

Fully numerical results for black-hole head-on collisions 
obtained in the equal and comparable mass regime have 
been compared with PP predictions and results obtained 
in the close-limit approximation [s^] by Anninos and col- 
laborators [si, HI] . These studies demonstrated agree- 
ment for the radiated energy and linear momentum bear- 
ing in mind the accuracies achievable at the time. How- 
ever, the spurious radiation present in the initial data was 
not dealt with in a satisfactory manner. Its presence con- 
taminated the physical pulse and much of the conclusions 
in these earlier works were affected by this. In particular, 
the presence of spurious radiation prevented an accurate 
computation of the total radiation and comparison with 
the linearized, PP calculations. 

Also, we are not aware of any comparisons between PP 
calculations and fully numerical results for mass ratios in 
a truly perturbative regime. By simulating black-hole bi- 
naries up to a mass ratio oi q = 1/100 we fill this gap 
and identify those aspects of the PP predictions which 
describe black-hole dynamics well in general and which 
only hold in the extreme mass-ratio limit. From a differ- 
ent point of view, the agreement with the PP calculations 
represents an important validation of the fully numerical 
calculations in the regime of high-mass ratios. In this 
context we emphasize that we are able to accurately ex- 
tract from binary black-hole simulations radiated GW 
energies of the order of 10~^ M and linear momenta cor- 
responding to recoil velocities of a few dozens of m/s, 
similar to the average speed of a normal car. We note, 
however, that even smaller amounts of energy have been 
extracted from general relativistic simulations of stellar 
core collapse; see e. g. [ssl ]. 

This paper is organized as follows. We summarize our 
numerical framework in Sec. II, estimate numerical un- 
certainties in Sec. Ill, describe our results in Sec. IV and 
conclude in Sec. V. 



II. NUMERICAL SETUP AND ANALYSIS 
TOOLS. 

The numerical simulations of unequal-mass black-hole 
collisions starting from rest have been performed with the 
Lean code, originally introduced in Ref. [s^, [s^]. The 
Lean code is based on the Cactus computational toolkit 
381 l39l and uses the Carpet mesh refinement package 
tl4l|. the apparent horizon finder AHFinderDirect [42|, 
l43j | and the TwoPuncture initial data solver [44|. The 
3-1-1 Einstein's equations are evolved using the BSSN 
[isl . formulation, together with the moving puncture 
approach 0, Q • The gauge conditions are determined by 
the so-called puncture gauge, i.e., the "1-l-log" slicing and 
r driver shift condition |47| . The systems are set up using 
Brill-Lindquist initial data. We have evolved BH binaries 
with mass ratios q = TO2/mi = 1, 1/2, 1/3, 1/4, 1/10 and 
1/100, where rrii is the bare mass parameter of the i-th 
BH. 

We use the Newman-Penrosc scalar ^'4 to measure 
gravitational radiation at extraction radii Rex, chosen in 
a range of 40 M to 90 M from the center of the colli- 
sion. We decompose ^'4 into multipoles ipim using spher- 
ical harmonics of spin- weight —2, _2llm, according to 

Due to the symmetry properties of the systems under 
consideration, the only non-vanishing multipoles all have 
TO = in a suitably chosen frame, and are purely real, 
corresponding to a single polarization state . In the 
equal-mass limit, the additional symmetry causes all mul- 
tipoles with odd I to vanish identically. The energy spec- 
trum and luminosity of the radiation are given by 
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respectively, where a hat denotes the Fourier transform 
and ipio is evaluated on a sphere at infinity. 



III. SIMULATIONS AND UNCERTAINTIES 

We have performed a series of simulations of head- 
on collisions with mass ratio ranging from q = 1 to 
q = 1/100 with initial coordinate separation d and proper 
horizon-to-horizon separation L as given in Table HI We 
describe the grid setup used for these simulations in terms 
of the number riii of refinement levels, the radius R of 
the computational domain, the resolution H used in the 
wave extraction zoncQ, the radius r in units of the smaller 
hole's mass TO2 of the innermost refinement level centered 



^ Typically the third refinement level counted from the outside. 
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TABLE I. Mass ratio q, coordinate and proper separation d 
and L, respectively, as well as radiated energy iSrad with per- 
centage distribution in the Z = 2, / = 3 and I = 4 muhipoles 
and recoil velocity v for the set of binary models evolved nu- 
merically. 
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TABLE II. Grid setup used for the different mass ratios q. 
The number of refinement levels is given by n^i , R is the radius 
of the computational domain, H the resolution in the wave 
extraction zone, r the radius of the innermost refinement box 
around the individual punctures and h the resolution used on 
that level. The additional low and high resolution for q = 1/4 
and q — 1/100 have been used for the convergence studies. 



on the individual puncturetH and the resolution hjmi of 
the innermost refinement level. The values for these pa- 
rameters are summarized for all mass ratios in Table HIl 

Our results are affected by three main sources of un- 
certainties: finite extraction radius, discretization and, 
for small initial separations of the binary, spurious initial 
radiation. We reduce the error arising from finite ex- 
traction radius by measuring the waveform components 
at several radii, and fitting them to an expression of the 
form ipim{r,t) = ■>pl^{t) +^\^{t)/r. The waveform "at 



^ For the small mass ratios q = 1/10 (1/100), the two (five) highest 
resolution boxes are placed around the small hole only to reduce 
computational cost. 
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FIG. 1. Convergence analysis for the I — 2 multipole of 
the gravitational wave signal for simulation g = 1/4, D = 
16.53 M (upper) and simulation q = 1/100, D = 9.58 (lower 
panel). In both cases we show the higher resolution differ- 
ences (solid black) together with the lower resolution result 
rescaled for second (dashed red lines) and fourth-order con- 
vergence (dotted blue lines). 



infinity" ^/;|^' {t) is the quantity reported throughout this 
work and used to calculate related quantities, such as the 
radiated energy. The uncertainty in this extrapolated 
value is estimated by performing a second fit including 

also a quadratic term ip'l^/i''^, and taking the difference 
between the first- and second-order fits. The resulting 
uncertainty increases as we decrease the mass ratio q and 
is 1 — 4 % for the total radiated energy and the I ~ 2 
waveform and energy, and 3 — 5 % for the subdominant 
multipoles and the radiated linear momentum. 

In order to estimate the discretization error of our sim- 
ulations, we have performed a convergence analysis for 
models (g = 1/4, L = 16.53 M) and {q = 1/100, L = 
9.58 M) using the three resolutions listed for these mass 
ratios in Table |lll The resulting convergence plots for the 
I — 2 multipole of the wave signal is shown in Fig. [T]and 
demonstrates convergence between second and fourth or- 
der. With regard to the analysis below, we note in partic- 
ular that the q = 1/100 case exhibits second order conver- 
gence in the plunge- merger signal around t — i?ex ~ 40 M 
but is close to fourth-order convergence for the remainder 
of the waveform. Bearing in mind that the plunge- merger 
transition represents the most dynamic part of the evo- 
lution and that the second-order ingredients in the code 
are associated with the prolongation of grid functions 
at the refinement boundaries in time, this observation 
is compatible with the numerical discretization. We ob- 
serve similar convergence properties for the I = 3 mul- 
tipole, but overall convergence close to fourth-order for 
the radiated energy and linear momentum, presumably 
because the accumulated errors are dominated by the 
fourth-order contributions observed for most of the sig- 
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nal. The resulting numerical uncertainties for q = 1/100 
are about 10 % in the waveform for the plunge-merger 
transition and 5 % for the remainder of the signal as well 
as 6 % for the radiated energy and 8 % for the linear 
momentum lost in gravitational waves. We note that in 
both cases, the discretization error leads to an overesti- 
mate of the radiated quantities. For g = 1/4 we observe 
significantly smaller uncertainties in the range of 2 % for 
all quantities. 

Finally, we comment on the unphysical gravitational 
radiation inherent in the conformally flat puncture ini- 
tial data. In order to extract physically meaningful infor- 
mation, one has to separate the spurious radiation from 
the radiation generated by the collision itself. This is 
done by "waiting" for the spurious radiation to radiate 
off the computational domain, and then discarding the 
early, contaminated part of the wave signal. For small 
values of the initial separation, however, the binary will 
merge before the spurious radiation has had enough time 
to leave the system, and physical and unphysical contri- 
butions to the wave signal partially overlap and cannot 
be cleanly distinguished. For our set of simulations, this 
problem arises only in the case q = 1/100, L ~ 9.58 M, 
where it introduces an additional error of about 2 % to 
the radiated energy and momentum. 



IV. RESULTS 

All collisions summarized in TableUresult in the forma- 
tion of a single BH plus gravitational radiation, i. e. there 
is no indication of violation of the cosmic censorship con- 
jecture. The final BH is born distorted, and eventually 
rings down to a Schwarzschild solution via emission of a 
superposition of quasi- normal modes [48 1. 

We illustrate the / = 2 and / = 3 wave signal in Fig. [2] 
for the / ~ 2 and I ~ 3 multipoles obtained for the mass 
ratios q = 1/4 (top), q = 1/10 (center) and q = 1/100 



(bottom) . In each panel the solid (black) curves represent 
the PP prediction for infall from infinity whereas the dot- 
ted (red) and dash-dotted (blue) curves show the numer- 
ical results for different values of the finite initial separa- 
tion. To leading order, the gravitational radiation output 
of black-hole collisions scales with the square of the re- 
duced mass fj. = Mrj of the system, where rj = (?/(<7 + 1)^ 
is the dimensionless, symmetric mass ratio [23|. For com- 
parison of the numerical results with PP predictions, we 
therefore rescale the former by the corresponding powers 
of 77, quadratic for energy and linear for the waveforms 
in Fig. H 

The waveforms show interesting features. For small 
initial separations, the early part of the waveform is con- 
taminated by "spurious" radiation; cf. the dotted (red) 
curve in the top and bottom panels of Fig. [5] As the 
initial separation increases, however, this problem disap- 
pears, because the longer infall duration of the binary 
provides sufficient time for the unphysical radiation to 
propagate off the grid; cf. the dash-dotted (blue) curves. 
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FIG. 2. (Color online) Rescaled waveforms for mass ratios q = 
1/4 (top), q = 1/10 (center) and q — 1/100 (bottom panels) 
for / = 2 (upper) and I = 3 (lower half of each panel), for two 
different initial separations. Also shown is the waveform in 
the PP limit (black solid lines). 
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vL (km/s) 



1/1 1/4 1/10 1/100 PP 
0.00936 0.00911 0.00985 0.0114 0.0104 
0.0 258.0 250.3 275.9 257.6 



TABLE III. Summary of our results when fitted to Eqs. (|4]) 
and The last column refers to PP results, as extrapolated 
from Lousto and Price [281. 



A closer inspection of the q = 1/100 case yields excellent 
agreement between the numerical and PP predictions ex- 
cept for the plunge-merger transition around < « in the 
figure. From the discussion in Sec. Ill, however, we re- 
call that the discretization error is particularly large in 
this regime. In fact, for the q = 1/100 model studied 
in Sec. Ill, a second-order Richardson extrapolation pre- 
dicts about a 10 % reduction in the amplitude around 
the first strong maximum in the I = 2 waveform which 
is very close in magnitude and sign to the deviation of 
the numerical from the PP result. As demonstrated by 
the upper central panel in Fig. [2l we find equally good 
agreement of the numerical / = 2 multipole with PP pre- 
dictions for the less extreme mass ratio g = 1/10 and only 
a small deviation for the larger mass ratio g = 1/4 (up- 
per top panel in Fig. [2]). Our findings thus confirm over 
a wide range of mass ratios the observation by Ref. [11] , 
that there is a weak dependence of the re-scaled wave- 
forms on the mass ratio. The I = 3 mode, on the other 
hand, is a good discriminator between high- and low- 
mass ratios. This behavior was qualitatively expected, 
as higher multipoles are suppressed in the equal-mass 
case; by symmetry the I — 3 mode is absent when the 
masses are equal. It is interesting, however, that even for 
what one might call a small mass ratio, q ~ 1/10, higher 
multipoles are still visibly suppressed. 

The total amount of energy radiated in gravitational 
waves during the collision depends on the initial separa- 
tion of the holes. As discussed in Anninos et al. [s^l, two 
effects contribute to increasing the GW energy at larger 
initial separations; (i) there is more time to radiate GWs 
during the infall and (ii) the infalling velocity is larger. 
In practice, the second effect is found to be dominant. 
Anninos et al. have accounted for both contributions by 
defining 
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One can write the corrections to the radiation emission 
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With the above as motivation, we have fitted our results 
to a 1/L dependence, of the form 
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with E^'^ the radiated energy for infinite initial sepa- 
ration. The results are summarized in Table IIIII We 
remind the reader that L stands for proper initial sepa- 
ration between the holes. We also note that the results 
in Table IIIII are normalized by rj^ . For comparison, we 
also show in the last entry of the table the results ob- 
tained in the PP limit, within a linearized calculation. 
This study was done by Lousto and Price [l^ using the 
same type of initial data; we have used their Table I to 
obtain the behavior shown in Table IIIII above. We note 
that already for q = 1/10 and q = 1/100 our results are 
in good agreement with PP calculations. We remind the 
reader, however, that in the q = 1/10 case there is a 
larger deviation in the / = 3 modes. 

With the extrapolation above one gets an estimate for 
the total radiation of two black holes merging from in- 
finite initial separation. A best fit of this number as 
function of mass ratio yields 



rad 



= 0.0110 - 0.00 8877 



(5) 



In the PP limit, when 77 —J- 0, this agrees with the clas- 
sical PP calculation, Eq. ([1]) to within 6%, so within the 
numerical uncertainties. Overall, the results in Table 
|l] demonstrate that we are able to accurately measure 
amounts of order E'^^'^ ~ 10~^M in these fully nonlinear 
evolutions. 

The amount of spurious radiation in the initial data 
is also consistent with predictions from linearized grav- 
ity. Lousto and Price performed a detailed analysis of 
the amount of spurious radiation in the infall of PPs into 
massive black holes, using the same type of initial data 
[28t . Using their Table I for L > 11, we find that the 
amount of spurious radiation varies with L according to 
S,.ad/(A/?72) - 0.15(L/M)-2 -5. For q = 1/100, for in- 
stance, we obtain E,^dl{Mrf) = 0.26(L/M)-2-55. Thus, 
we find good agreement in the decay power (roughly 
—2.5) and also in the proportionality coefficient. 

If two BHs with different masses collide head-on, the 
remnant BH will recoil with respect to the center-of-mass 
frame, due to the emission of energy and momentum car- 
ried by gravitational waves. Based on PN tools, we have 
fit our results to [i^ 

^^rocoil = (1 + M/L) , (6) 

where is a normalized recoil velocity for infinite initial 
separation. The normalized recoil velocity is shown 
in Table IIIII The point particle limit was considered in 
Ref. [iOl, who obtained = 263km/sE|. We note this 
is not a trivial agreement: unlike energy calculations, 
momentum involves interference with higher (typically 



^ note the sUght disagreement with the extrapolation of Lousto 
and Price's results, shown in Table UlTl 
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highly suppressed) multipoles. Overall, our results agree 
well in the limit of small mass-ratios with the point par- 
ticle limit. It is interesting to note in this context that 
for both, radiated energy and linear momentum, the nu- 
merical results exceed those obtained from the point par- 
ticle limit by about 6 %. This value agrees in sign and 
magnitude with the discretization error obtained for the 
q ~ 1/100 simulation in Sec. IIIII We therefore consider 
the discretization error the dominant source of the re- 
maining discrepancies. 



V. CONCLUSIONS 

The simulation of dynamical, interacting black holes 
has a tremendous potential to provide answers to some 
of the most fundamental questions in physics. Recent de- 
velopments in experimental and theoretical physics make 
this a pressing issue. We refer, in particular, to the 
prominent role of BHs in the gauge-gravity duality, in 
TeV-scale gravity or even on their own as solutions of 
the field equations [5l| . Recent work along these lines in- 
cludes the successful simulation and understanding of the 
collision of two BHs at close to the speed of light in four- 
dimensional spacetime [52| - [55| . the low energy collisions 
in higher spacetime dimensions [5l|,[5^,[53l, BH scattering 
in five dimensions [58| , stability studies in higher dimen- 
sions [59l - [6l| and BH evolutions in non asymptotically 
flat spacetimes (62| (for the formalism extension, we re- 
fer the reader to Refs. HH, M, M^)- 

We have shown here that Numerical Relativity is ca- 
pable of simulating dynamical black holes close to the 
regime of validity of linear calculations, and to make 
contact with approximation techniques. For this purpose 
we have evolved head-on collisions of non-spinning black- 
hole binaries over a range of mass ratios from q = 1 to 
q = 1/100. We obtain radiated energies decreasing from 
about 5.5 X 10"'^ for g = 1 to 10"^ for q = 1/100. The 
recoil reaches a maximum of about 4 km/s near q = 3 
and decreases towards 26 m/s for q = 1/100. In the 
limit of small mass ratios and extrapolating our results 
to infinite initial separation, we find the numerical values 
for radiated energy and linear momentum to be « 6 % 
larger than the point-particle predictions. This discrep- 
ancy agrees rather well in sign and magnitude with the 
discretization error obtained from a convergence study 



of our q = 1/100 simulations. It thus appears likely 
that a significant part of the remaining differences can 
be attributed to the discretization error which mirrors 
the computational demands of numerical black-hole bi- 
nary simulations with such small mass ratios. 

With regard to the waveforms, the most remarkable 
result is the suppression of odd / multipoles. While we 
observe good agreement between numerical and point- 
particle results for the / = 2 mode, already for q = 1/10, 
the numerically calculated I = 3 multipole is visibly sup- 
pressed for this case and only agrees well with the PP 
limit for q = 1/100. 

Overall, the good agreement for waveforms and radi- 
ated energy and momenta for the case q = 1/100 demon- 
strates that numerical techniques are capable of bridging 
the gap between linear analysis and the fully non-linear 
regime of general relativity. 
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